library(tidyverse)
library(patchwork)
b2018 <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/b2018.csv", header = TRUE)
b2019 <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/b2019.csv", header = TRUE)
b2020 <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/b2020.csv", header = TRUE)
bg <- read.csv("C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/datos/base/Chiapas 5/summer/data/nd/bg.csv", header = TRUE)Ahora creemos 3 graficos para cada año correspondiente de nuestro conjunto de datos, cargemos algunas librerias necesarias
La siguiente forma de mapear nuestro conjunto de datos obtenido, es de la siguiente forma, considere que esta forma no involucra leer un shape file, convertir los datos a tipo geodata y todas esas particularidades, eso es lo interesante y practico, aunque no descartamos hacerlo de la forma correspondiente.
Algo particular es que si necesitamos el shape file, solo para colocar el contorno del estado de chiapas, es por ello que debemos de leerlo.
cargamos algunas librerias
library(rgdal)
library(dplyr)
library(ggplot2)
library(leaflet) # libreria para graficar mapas interactivos
library(sf) # manejo de informacion geografica
library(viridis) # paletas de colores
library(RColorBrewer) # mas paletas de colores
library(patchwork)a continuación leemos nuestro archivo shapefile, algo que debemos identificar del shapefile es que su proyeccion geografica no coincide con la proyeccion geografica del conjunto de daos a trabajar, es por ello que se necesita realizar una transformación en la proyección de nuestro shapefile, para llevarlo a la proyección habitual de longitud y latitud que se conoce.
my_spdf <- readOGR(
dsn= "C:/Users/David/OneDrive/Documentos/MMA/seminario de tesis 2/avances seminario/A2",
layer="ENTIDAD",
verbose=FALSE
)
# trasformamos a el sistema de coordenadas habitual
my_spdf <- spTransform(my_spdf, CRS("+proj=longlat +datum=WGS84"))y seleccionamos el estado de chiapas
my_spdf_c <- my_spdf[my_spdf$nombre == "CHIAPAS",]podemos observar mediante un grafico que el estado seleccionado sea el correcto.
plot(my_spdf_c, col="#f2f2f2", bg="skyblue", lwd=0.25)Ahora si estamos listos para realizar un grafico de nuestro conjunto de datos
b2018 %>%
ggplot() +
geom_point(aes(x = Longitud, y = Latitud, colour = Rain),size =3)+
borders(my_spdf_c)+
coord_quickmap()+
theme_test()Agregando un poco mas de detalle podemos obtener un grafico mas detallado para cada año y general
library(ggrepel)
b2018 %>%
ggplot(aes(x = Longitud, y = Latitud)) +
borders(my_spdf_c, fill = "antiquewhite1") +
geom_point(aes(colour = Rain), size = 4) +
scale_color_viridis_c(option = "plasma", trans = "sqrt",
oob = scales::squish) +
coord_quickmap() +
theme_test() +
xlab("Longitude") + ylab("Latitude") +
ggtitle("Chiapas 2018", subtitle = "Precipitación media por hora") +
theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed",
size = 0.5), panel.background = element_rect(fill = "aliceblue")) +
geom_text_repel(data = distinct(b2018, Longitud, .keep_all = TRUE),
aes(x = Longitud, y = Latitud, label = Station),
box.padding = 1, point.padding = 1, segment.color = 'grey50',
size = 4, point.size = 4, segment.size = 1) -> p5;p5# unique(b2018$Station)
# distinct(b2018, Longitud, .keep_all = T)
# table(b2018$Latitud)
b2019 %>%
ggplot() +
borders(my_spdf_c, fill = "antiquewhite1") +
geom_point(aes(x = Longitud, y = Latitud, color = Rain, label=Station), size =4) +
geom_text_repel(data = distinct(b2019, Station, .keep_all = TRUE),
aes(x = Longitud, y = Latitud, label = Station),
box.padding = 1, point.padding = 1, segment.color = 'grey50',
fontface = "bold", size = 4, point.size = 4, segment.size = 1)+
coord_quickmap() +
theme_test() +
scale_color_viridis_c(option = "plasma", trans = "sqrt") +
xlab("Longitude") + ylab("Latitude") +
ggtitle("Chiapas 2019", subtitle = "Precipitación media por hora") +
theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", size = 0.5),
panel.background = element_rect(fill = "aliceblue")) -> p6;p6b2020 %>%
ggplot() +
borders(my_spdf_c, fill = "antiquewhite1") +
geom_point(aes(x = Longitud, y = Latitud, color = Rain, label = Station), size =4) +
geom_text_repel(data = distinct(b2020, Station, .keep_all = TRUE),
aes(x = Longitud, y = Latitud, label = Station),
box.padding = 1, point.padding = 1, segment.color = 'grey50',
fontface = "bold", size = 4, point.size = 4, segment.size = 1)+
coord_quickmap() +
theme_test() +
scale_color_viridis_c(option = "plasma", trans = "sqrt") +
xlab("Longitude") + ylab("Latitude") +
ggtitle("Chiapas 2020", subtitle = "Precipitación media por hora") +
theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed",
size = 0.5), panel.background = element_rect(fill = "aliceblue")) -> p7;p7bg %>%
ggplot() +
borders(my_spdf_c, fill = "antiquewhite1") +
geom_point(aes(x = Longitud, y = Latitud, color = Rain, label=Station), size =4) +
geom_text_repel(data = distinct(b2020, Station, .keep_all = TRUE),
aes(x = Longitud, y = Latitud, label = Station),
box.padding = 1, point.padding = 1, segment.color = 'grey50',
fontface = "bold", size = 4, point.size = 4, segment.size = 1)+
coord_quickmap() +
theme_test() +
scale_color_viridis_c(option = "plasma", trans = "sqrt") +
xlab("Longitude") + ylab("Latitude") +
ggtitle("Chiapas 2018-2020", subtitle = "Precipitación media por hora") +
theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed",
size = 0.5), panel.background = element_rect(fill = "aliceblue")) -> p8;p8(p5 + p6)/(p7 + p8)Para este caso, lo que necesitamos son ubicaciones en el estado de Chiapas, esto con el objetivo de poder simular en esas ubicaciones el modelo de valores extremos, para ello vamos utilizar la libreria geoR, esto con el fin de predecir en este caso 40 ubicaciones en la zona de estudio.
library(geoR)
bor1 <- my_spdf_c@polygons[[1]]@Polygons[[1]]@coords
sim1 <- grf(n= 40, grid = "irreg", cov.pars=c(2, 1), borders = bor1, nugget = 5, mean = 200,
cov.model = "gaussian")> grf: process with 1 covariance structure(s)
> grf: nugget effect is: tausq= 5
> grf: covariance model 1 is: gaussian(sigmasq=2, phi=1)
> grf: decomposition algorithm used is: cholesky
> grf: End of simulation procedure. Number of realizations: 1
Podemos observar el mapa de las ubicaciones de interes simuladas, estas ubicaciones se encuentran aleatoriamente en la zona de estudio, lo unico que nos interesa en este caso es la ubicación de estos puntos simulados
par(mar=c(3,3,3,0))
points(sim1, main ="Chiapas SIM")Conforme a la metodologia de Davison, A.C., Padoan, S.A. and Ribatet. Lo que realizo es en los 40 sitios generados aleatoriamente, voy a obtener 50 observaciones en cada uno de ellos mediante la generación de valores de una
library(SpatialExtremes)
## Not run:
## Generate realizations from the model
n.site <- 40
n.obs <- 50
# class(sim1$coords)
coord <- cbind(lon = sim1$coords[,1], lat = sim1$coords[,2])
gp.loc <- rgp(1, coord, "powexp", sill = 4, range = 20, smooth = 1)
gp.scale <- rgp(1, coord, "powexp", sill = 4, range = 20, smooth = 1)
gp.shape <- rgp(1, coord, "powexp", sill = 2, range = 20, smooth = 1)
# locs <- 26 + 0.5 * coord[,"lon"] + gp.loc
# scales <- 10 + 0.2 * coord[,"lat"] + gp.scale
# shapes <- 0.15 + gp.shape
locs <- gp.loc
scales <- abs(gp.scale)
shapes <- gp.shape
data <- matrix(NA, n.obs, n.site)
for (i in 1:n.site)
data[,i] <- rgev(n.obs, locs[i], scales[i], shapes[i])
loc.form <- y ~ 1
scale.form <- y ~ 1
shape.form <- y ~ 1
hyper <- list()
hyper$sills <- list(loc = c(1,8), scale = c(1,1), shape = c(1,0.02))
hyper$ranges <- list(loc = c(2,20), scale = c(1,5), shape = c(1, 10))
hyper$smooths <- list(loc = c(1,1/3), scale = c(1,1/3), shape = c(1, 1/3))
hyper$betaMeans <- list(loc = 20,
scale = 15,
shape = 10)
hyper$betaIcov <- list(loc = solve(diag(c(10), 1, 1)),
scale = solve(diag(c(10), 1, 1)),
shape = solve(diag(c(10), 1, 1)))
## We will use an exponential covariance function so the jump sizes for
## the shape parameter of the covariance function are null.
prop <- list(gev = c(1.2, 0.08, 0.08), ranges = c(0.7, 0.8, 0.7), smooths = c(0,0,0))
start <- list(sills = c(4, .36, 0.009),
ranges = c(24, 17, 16),
smooths = c(1, 1, 1),
beta = list(loc = c(0.2),
scale = c(0.3),
shape = c(20)))
# mc <- latent(data, coord, loc.form = loc.form, scale.form = scale.form,
# shape.form = shape.form, hyper = hyper, prop = prop, start = start,
# n = 10000, burn.in = 5000, thin = 15)
# mc
mc1 <- latent(data, coord,
loc.form = loc.form,
scale.form = scale.form,
shape.form = shape.form,
hyper = hyper,
prop = prop,
start = start,
n = 10000,
burn.in = 5000,
thin = 15)
mc1> Effective length: 10000
> Burn-in: 5000
> Thinning: 15
> Effective NoP: 57.72025
> DIC: 10420.48
>
> Regression Parameters:
> Location Parameters:
> lm1
> ci.lower 7.986
> post.mean 14.230
> ci.upper 20.550
>
> Scale Parameters:
> lm1
> ci.lower 1.604
> post.mean 3.519
> ci.upper 7.058
>
> Shape Parameters:
> lm1
> ci.lower 0.2127
> post.mean 0.7418
> ci.upper 1.8975
>
>
> Latent Parameters:
> Location Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 8.346 33.328 1.000
> post.mean 25.614 90.640 1.000
> ci.upper 58.684 178.641 1.000
>
> Scale Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 0.3121 2.3281 1.0000
> post.mean 1.4136 10.4406 1.0000
> ci.upper 4.8123 26.8158 1.0000
>
> Shape Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 0.02498 0.25977 1.00000
> post.mean 0.19374 3.59892 1.00000
> ci.upper 1.06934 19.52523 1.00000
mc2 <- latent(data, coord,
loc.form = loc.form,
scale.form = scale.form,
shape.form = shape.form,
hyper = hyper,
prop = prop,
start = start,
n = 10000,
burn.in = 5000,
thin = 15)
mc2> Effective length: 10000
> Burn-in: 5000
> Thinning: 15
> Effective NoP: 57.72025
> DIC: 10420.48
>
> Regression Parameters:
> Location Parameters:
> lm1
> ci.lower 7.986
> post.mean 14.230
> ci.upper 20.550
>
> Scale Parameters:
> lm1
> ci.lower 1.604
> post.mean 3.519
> ci.upper 7.058
>
> Shape Parameters:
> lm1
> ci.lower 0.2127
> post.mean 0.7418
> ci.upper 1.8975
>
>
> Latent Parameters:
> Location Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 8.346 33.328 1.000
> post.mean 25.614 90.640 1.000
> ci.upper 58.684 178.641 1.000
>
> Scale Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 0.3121 2.3281 1.0000
> post.mean 1.4136 10.4406 1.0000
> ci.upper 4.8123 26.8158 1.0000
>
> Shape Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 0.02498 0.25977 1.00000
> post.mean 0.19374 3.59892 1.00000
> ci.upper 1.06934 19.52523 1.00000
mc3 <- latent(data, coord,
loc.form = loc.form,
scale.form = scale.form,
shape.form = shape.form,
hyper = hyper,
prop = prop,
start = start,
n = 10000,
burn.in = 5000,
thin = 15)
mc3> Effective length: 10000
> Burn-in: 5000
> Thinning: 15
> Effective NoP: 57.72025
> DIC: 10420.48
>
> Regression Parameters:
> Location Parameters:
> lm1
> ci.lower 7.986
> post.mean 14.230
> ci.upper 20.550
>
> Scale Parameters:
> lm1
> ci.lower 1.604
> post.mean 3.519
> ci.upper 7.058
>
> Shape Parameters:
> lm1
> ci.lower 0.2127
> post.mean 0.7418
> ci.upper 1.8975
>
>
> Latent Parameters:
> Location Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 8.346 33.328 1.000
> post.mean 25.614 90.640 1.000
> ci.upper 58.684 178.641 1.000
>
> Scale Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 0.3121 2.3281 1.0000
> post.mean 1.4136 10.4406 1.0000
> ci.upper 4.8123 26.8158 1.0000
>
> Shape Parameters:
> Covariance family: powexp
> sill range smooth
> ci.lower 0.02498 0.25977 1.00000
> post.mean 0.19374 3.59892 1.00000
> ci.upper 1.06934 19.52523 1.00000
dsgn.mat The design matrix.
sill = umbral
ranges = rango
smooths = forma
Un vector de longitud 3 que devuelve el DIC, el número efectivo de parámetros eNoP y una estimación de la desviación esperada Dbar.
library(broom)
library(coda)
library(broom.mixed)
library(brms)
library(bayesplot)
library(lattice)library(lattice)
mc1$chain.loc <- as.data.frame(mc1$chain.loc) %>%
mutate(chain = 1)
mc2$chain.loc <- as.data.frame(mc2$chain.loc) %>%
mutate(chain = 2)
mc3$chain.loc <- as.data.frame(mc3$chain.loc) %>%
mutate(chain = 3)
mc.loc <- mc1$chain.loc %>%
union_all(mc2$chain.loc) %>% union_all(mc3$chain.loc)
post <- posterior_samples(mc.loc, add_chain = T)
mcmc_dens_overlay(post)post %>%
select(lm1, sill, range, loc1, loc2, loc3, loc4) %>%
mcmc_acf()post %>%
select(lm1, sill, range, loc1, loc2, loc3, loc4) %>%
mcmc_trace()post %>%
select(lm1, sill, range, loc1, loc2, loc3, loc4) %>%
mcmc_dens()post %>% select(lm1) %>% mcmc_trace()library(coda)
coda::raftery.diag(posterior_samples(post))>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3994 3746 1.07
> sill 3 4096 3746 1.09
> range 4 4673 3746 1.25
> smooth <NA> <NA> 3746 NA
> loc1 8 10268 3746 2.74
> loc2 4 5080 3746 1.36
> loc3 8 9798 3746 2.62
> loc4 6 8112 3746 2.17
> loc5 9 9495 3746 2.53
> loc6 4 4752 3746 1.27
> loc7 3 4061 3746 1.08
> loc8 5 6024 3746 1.61
> loc9 3 4373 3746 1.17
> loc10 4 4913 3746 1.31
> loc11 9 12591 3746 3.36
> loc12 4 4884 3746 1.30
> loc13 4 4831 3746 1.29
> loc14 4 4831 3746 1.29
> loc15 4 5254 3746 1.40
> loc16 5 5390 3746 1.44
> loc17 3 4061 3746 1.08
> loc18 4 5038 3746 1.34
> loc19 6 6462 3746 1.73
> loc20 10 12176 3746 3.25
> loc21 3 4232 3746 1.13
> loc22 4 4996 3746 1.33
> loc23 8 8572 3746 2.29
> loc24 6 8466 3746 2.26
> loc25 6 7200 3746 1.92
> loc26 6 8462 3746 2.26
> loc27 8 9006 3746 2.40
> loc28 7 7400 3746 1.98
> loc29 4 4673 3746 1.25
> loc30 4 4673 3746 1.25
> loc31 6 8036 3746 2.15
> loc32 6 7656 3746 2.04
> loc33 8 11482 3746 3.07
> loc34 6 9066 3746 2.42
> loc35 15 17613 3746 4.70
> loc36 6 8122 3746 2.17
> loc37 4 4597 3746 1.23
> loc38 4 4831 3746 1.29
> loc39 4 4752 3746 1.27
> loc40 8 9117 3746 2.43
> chain 69075 69075 3746 18.40
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.loc[,c(-4,-9)]),
# as.mcmc(mc2$chain.loc[,c(-4,-9)]),
# as.mcmc(mc3$chain.loc[,c(-4,-9)])))
#
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.loc[,1]),
# as.mcmc(mc2$chain.loc[,1]),
# as.mcmc(mc3$chain.loc[,1])))
BMu1.mcmc<-mcmc.list(as.mcmc(mc1$chain.loc[,c(-4,-9)]),
as.mcmc(mc2$chain.loc[,c(-4,-9)]),
as.mcmc(mc3$chain.loc[,c(-4,-9)]))
summary(BMu1.mcmc)>
> Iterations = 1:10000
> Thinning interval = 1
> Number of chains = 3
> Sample size per chain = 10000
>
> 1. Empirical mean and standard deviation for each variable,
> plus standard error of the mean:
>
> Mean SD Naive SE Time-series SE
> lm1 14.2297 3.2277 0.0186353 0.020648
> sill 25.6144 12.8728 0.0743212 0.093081
> range 90.6396 37.5718 0.2169210 0.268302
> loc1 1.6649 0.1669 0.0009634 0.001768
> loc2 1.1084 0.2365 0.0013654 0.002141
> loc3 2.6815 0.1991 0.0011494 0.002124
> loc4 1.4826 0.1550 0.0008947 0.001900
> loc6 1.5297 0.2102 0.0012134 0.002090
> loc7 2.3124 0.3217 0.0018574 0.002499
> loc8 1.1358 0.2054 0.0011857 0.001996
> loc9 2.2184 0.1945 0.0011229 0.001706
> loc10 2.0980 0.1892 0.0010926 0.001660
> loc11 1.8553 0.1551 0.0008955 0.001669
> loc12 2.8993 0.2871 0.0016574 0.002579
> loc13 1.3152 0.1928 0.0011133 0.001533
> loc14 2.1239 0.2245 0.0012962 0.001910
> loc15 0.9606 0.1953 0.0011275 0.001722
> loc16 0.8013 0.2432 0.0014039 0.001933
> loc17 2.0278 0.1857 0.0010719 0.001463
> loc18 1.5992 0.2115 0.0012212 0.002068
> loc19 1.7372 0.1574 0.0009088 0.001827
> loc20 1.9681 0.1521 0.0008783 0.002631
> loc21 1.8076 0.2176 0.0012565 0.001625
> loc22 2.7312 0.2012 0.0011615 0.001811
> loc23 2.4389 0.2191 0.0012649 0.002205
> loc24 1.4109 0.2311 0.0013343 0.002245
> loc25 1.7159 0.1478 0.0008534 0.001889
> loc26 1.4967 0.2041 0.0011787 0.001733
> loc27 1.4447 0.1708 0.0009860 0.001962
> loc28 1.7872 0.1631 0.0009415 0.002100
> loc29 1.2776 0.3351 0.0019347 0.002832
> loc30 2.6040 0.2886 0.0016663 0.002571
> loc31 1.3298 0.1946 0.0011238 0.001751
> loc32 1.6352 0.1847 0.0010664 0.002063
> loc33 1.7319 0.1546 0.0008929 0.001864
> loc34 1.6563 0.1839 0.0010620 0.001982
> loc35 1.9140 0.1555 0.0008977 0.002602
> loc36 2.9123 0.1927 0.0011128 0.001833
> loc37 2.7145 0.2101 0.0012130 0.001996
> loc38 1.4030 0.1941 0.0011207 0.001710
> loc39 0.7753 0.2160 0.0012473 0.001728
> loc40 1.6935 0.1554 0.0008973 0.001975
> chain 2.0000 0.8165 0.0047141 0.000000
>
> 2. Quantiles for each variable:
>
> 2.5% 25% 50% 75% 97.5%
> lm1 7.9863 12.0235 14.1928 16.4447 20.550
> sill 8.3460 16.5638 22.9614 31.7648 58.684
> range 33.3278 63.4518 85.1777 111.7662 178.641
> loc1 1.3388 1.5532 1.6633 1.7754 2.005
> loc2 0.6341 0.9545 1.1121 1.2688 1.564
> loc3 2.3052 2.5479 2.6767 2.8095 3.082
> loc4 1.1857 1.3775 1.4792 1.5869 1.794
> loc6 1.1307 1.3868 1.5268 1.6691 1.955
> loc7 1.7000 2.0965 2.3065 2.5208 2.961
> loc8 0.7305 0.9979 1.1367 1.2741 1.535
> loc9 1.8506 2.0848 2.2127 2.3483 2.609
> loc10 1.7318 1.9717 2.0970 2.2243 2.480
> loc11 1.5622 1.7490 1.8529 1.9575 2.167
> loc12 2.3618 2.7045 2.8911 3.0802 3.503
> loc13 0.9393 1.1859 1.3147 1.4411 1.705
> loc14 1.6874 1.9732 2.1214 2.2720 2.573
> loc15 0.5750 0.8288 0.9608 1.0926 1.343
> loc16 0.3258 0.6391 0.7956 0.9628 1.285
> loc17 1.6709 1.9015 2.0237 2.1529 2.392
> loc18 1.1934 1.4570 1.5992 1.7395 2.024
> loc19 1.4373 1.6305 1.7339 1.8427 2.052
> loc20 1.6848 1.8640 1.9645 2.0684 2.279
> loc21 1.3817 1.6618 1.8072 1.9536 2.233
> loc22 2.3526 2.5925 2.7260 2.8648 3.140
> loc23 2.0164 2.2882 2.4370 2.5840 2.879
> loc24 0.9653 1.2586 1.4071 1.5629 1.865
> loc25 1.4339 1.6136 1.7136 1.8143 2.014
> loc26 1.1085 1.3571 1.4948 1.6317 1.913
> loc27 1.1114 1.3305 1.4411 1.5586 1.782
> loc28 1.4788 1.6756 1.7821 1.8958 2.125
> loc29 0.6158 1.0549 1.2772 1.5007 1.942
> loc30 2.0468 2.4090 2.5997 2.7892 3.180
> loc31 0.9537 1.1987 1.3286 1.4592 1.719
> loc32 1.2736 1.5121 1.6367 1.7608 2.004
> loc33 1.4307 1.6275 1.7295 1.8322 2.040
> loc34 1.2974 1.5318 1.6559 1.7798 2.015
> loc35 1.6182 1.8086 1.9090 2.0170 2.224
> loc36 2.5401 2.7812 2.9073 3.0403 3.304
> loc37 2.3147 2.5724 2.7103 2.8543 3.135
> loc38 1.0310 1.2713 1.4035 1.5297 1.791
> loc39 0.3605 0.6279 0.7756 0.9167 1.213
> loc40 1.3970 1.5887 1.6923 1.7956 2.007
> chain 1.0000 1.0000 2.0000 3.0000 3.000
xyplot(BMu1.mcmc)densityplot(BMu1.mcmc) #Densidadesplot(BMu1.mcmc)# layout(matrix(1:12, 3,4));
traceplot(BMu1.mcmc)#Grafica de autocorrelacion
autocorr.plot(BMu1.mcmc, auto.layout = TRUE, ask =F)coda::acfplot(BMu1.mcmc)# gelman.diag(BMu1.mcmc)
#
# gelman.plot(BMu1.mcmc)
geweke.diag(BMu1.mcmc)> [[1]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range loc1 loc2 loc3 loc4 loc6
> -2.02360 -0.13114 -0.57292 -0.13784 1.40118 -0.25933 0.79576 0.38023
> loc7 loc8 loc9 loc10 loc11 loc12 loc13 loc14
> -0.65418 -0.81671 0.76711 1.53365 -0.43782 -0.78151 0.43533 0.83240
> loc15 loc16 loc17 loc18 loc19 loc20 loc21 loc22
> 0.21704 -0.02015 2.09454 0.25178 -1.06965 -0.58544 1.15215 -2.70446
> loc23 loc24 loc25 loc26 loc27 loc28 loc29 loc30
> -1.37479 -0.25433 -0.39967 -0.52450 1.25959 -1.75285 -0.62833 -1.14631
> loc31 loc32 loc33 loc34 loc35 loc36 loc37 loc38
> -0.50914 0.96571 -0.59595 0.91762 -0.59011 0.74688 0.20454 -0.78983
> loc39 loc40 chain
> -1.10704 -0.34407 NaN
>
>
> [[2]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range loc1 loc2 loc3 loc4 loc6
> -2.02360 -0.13114 -0.57292 -0.13784 1.40118 -0.25933 0.79576 0.38023
> loc7 loc8 loc9 loc10 loc11 loc12 loc13 loc14
> -0.65418 -0.81671 0.76711 1.53365 -0.43782 -0.78151 0.43533 0.83240
> loc15 loc16 loc17 loc18 loc19 loc20 loc21 loc22
> 0.21704 -0.02015 2.09454 0.25178 -1.06965 -0.58544 1.15215 -2.70446
> loc23 loc24 loc25 loc26 loc27 loc28 loc29 loc30
> -1.37479 -0.25433 -0.39967 -0.52450 1.25959 -1.75285 -0.62833 -1.14631
> loc31 loc32 loc33 loc34 loc35 loc36 loc37 loc38
> -0.50914 0.96571 -0.59595 0.91762 -0.59011 0.74688 0.20454 -0.78983
> loc39 loc40 chain
> -1.10704 -0.34407 NaN
>
>
> [[3]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range loc1 loc2 loc3 loc4 loc6
> -2.02360 -0.13114 -0.57292 -0.13784 1.40118 -0.25933 0.79576 0.38023
> loc7 loc8 loc9 loc10 loc11 loc12 loc13 loc14
> -0.65418 -0.81671 0.76711 1.53365 -0.43782 -0.78151 0.43533 0.83240
> loc15 loc16 loc17 loc18 loc19 loc20 loc21 loc22
> 0.21704 -0.02015 2.09454 0.25178 -1.06965 -0.58544 1.15215 -2.70446
> loc23 loc24 loc25 loc26 loc27 loc28 loc29 loc30
> -1.37479 -0.25433 -0.39967 -0.52450 1.25959 -1.75285 -0.62833 -1.14631
> loc31 loc32 loc33 loc34 loc35 loc36 loc37 loc38
> -0.50914 0.96571 -0.59595 0.91762 -0.59011 0.74688 0.20454 -0.78983
> loc39 loc40 chain
> -1.10704 -0.34407 NaN
#geweke.plot(BMu1.mcmc)
raftery.diag(BMu1.mcmc)> [[1]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3994 3746 1.07
> sill 3 4098 3746 1.09
> range 4 4674 3746 1.25
> loc1 8 10270 3746 2.74
> loc2 4 5080 3746 1.36
> loc3 6 6349 3746 1.69
> loc4 6 6877 3746 1.84
> loc6 4 4752 3746 1.27
> loc7 3 4061 3746 1.08
> loc8 5 6025 3746 1.61
> loc9 3 4374 3746 1.17
> loc10 4 4913 3746 1.31
> loc11 4 4913 3746 1.31
> loc12 4 4884 3746 1.30
> loc13 4 4832 3746 1.29
> loc14 4 4832 3746 1.29
> loc15 4 5254 3746 1.40
> loc16 5 5390 3746 1.44
> loc17 3 4061 3746 1.08
> loc18 4 5038 3746 1.34
> loc19 6 6462 3746 1.73
> loc20 9 9691 3746 2.59
> loc21 3 4232 3746 1.13
> loc22 4 4996 3746 1.33
> loc23 5 6247 3746 1.67
> loc24 4 5166 3746 1.38
> loc25 6 7200 3746 1.92
> loc26 3 4520 3746 1.21
> loc27 6 6816 3746 1.82
> loc28 7 7401 3746 1.98
> loc29 4 4674 3746 1.25
> loc30 4 4674 3746 1.25
> loc31 4 5080 3746 1.36
> loc32 5 5820 3746 1.55
> loc33 7 7465 3746 1.99
> loc34 4 5254 3746 1.40
> loc35 10 11232 3746 3.00
> loc36 3 4520 3746 1.21
> loc37 4 4597 3746 1.23
> loc38 4 4832 3746 1.29
> loc39 4 4752 3746 1.27
> loc40 8 9118 3746 2.43
> chain <NA> <NA> 3746 NA
>
>
> [[2]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3994 3746 1.07
> sill 3 4098 3746 1.09
> range 4 4674 3746 1.25
> loc1 8 10270 3746 2.74
> loc2 4 5080 3746 1.36
> loc3 6 6349 3746 1.69
> loc4 6 6877 3746 1.84
> loc6 4 4752 3746 1.27
> loc7 3 4061 3746 1.08
> loc8 5 6025 3746 1.61
> loc9 3 4374 3746 1.17
> loc10 4 4913 3746 1.31
> loc11 4 4913 3746 1.31
> loc12 4 4884 3746 1.30
> loc13 4 4832 3746 1.29
> loc14 4 4832 3746 1.29
> loc15 4 5254 3746 1.40
> loc16 5 5390 3746 1.44
> loc17 3 4061 3746 1.08
> loc18 4 5038 3746 1.34
> loc19 6 6462 3746 1.73
> loc20 9 9691 3746 2.59
> loc21 3 4232 3746 1.13
> loc22 4 4996 3746 1.33
> loc23 5 6247 3746 1.67
> loc24 4 5166 3746 1.38
> loc25 6 7200 3746 1.92
> loc26 3 4520 3746 1.21
> loc27 6 6816 3746 1.82
> loc28 7 7401 3746 1.98
> loc29 4 4674 3746 1.25
> loc30 4 4674 3746 1.25
> loc31 4 5080 3746 1.36
> loc32 5 5820 3746 1.55
> loc33 7 7465 3746 1.99
> loc34 4 5254 3746 1.40
> loc35 10 11232 3746 3.00
> loc36 3 4520 3746 1.21
> loc37 4 4597 3746 1.23
> loc38 4 4832 3746 1.29
> loc39 4 4752 3746 1.27
> loc40 8 9118 3746 2.43
> chain <NA> <NA> 3746 NA
>
>
> [[3]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3994 3746 1.07
> sill 3 4098 3746 1.09
> range 4 4674 3746 1.25
> loc1 8 10270 3746 2.74
> loc2 4 5080 3746 1.36
> loc3 6 6349 3746 1.69
> loc4 6 6877 3746 1.84
> loc6 4 4752 3746 1.27
> loc7 3 4061 3746 1.08
> loc8 5 6025 3746 1.61
> loc9 3 4374 3746 1.17
> loc10 4 4913 3746 1.31
> loc11 4 4913 3746 1.31
> loc12 4 4884 3746 1.30
> loc13 4 4832 3746 1.29
> loc14 4 4832 3746 1.29
> loc15 4 5254 3746 1.40
> loc16 5 5390 3746 1.44
> loc17 3 4061 3746 1.08
> loc18 4 5038 3746 1.34
> loc19 6 6462 3746 1.73
> loc20 9 9691 3746 2.59
> loc21 3 4232 3746 1.13
> loc22 4 4996 3746 1.33
> loc23 5 6247 3746 1.67
> loc24 4 5166 3746 1.38
> loc25 6 7200 3746 1.92
> loc26 3 4520 3746 1.21
> loc27 6 6816 3746 1.82
> loc28 7 7401 3746 1.98
> loc29 4 4674 3746 1.25
> loc30 4 4674 3746 1.25
> loc31 4 5080 3746 1.36
> loc32 5 5820 3746 1.55
> loc33 7 7465 3746 1.99
> loc34 4 5254 3746 1.40
> loc35 10 11232 3746 3.00
> loc36 3 4520 3746 1.21
> loc37 4 4597 3746 1.23
> loc38 4 4832 3746 1.29
> loc39 4 4752 3746 1.27
> loc40 8 9118 3746 2.43
> chain <NA> <NA> 3746 NA
# heidel.diag(BMu1.mcmc)mc1$chain.scale <- as.data.frame(mc1$chain.scale) %>%
mutate(chain = 1)
mc2$chain.scale <- as.data.frame(mc2$chain.scale) %>%
mutate(chain = 2)
mc3$chain.scale <- as.data.frame(mc3$chain.scale) %>%
mutate(chain = 3)
mc.loc <- mc1$chain.scale %>%
union_all(mc2$chain.scale) %>% union_all(mc3$chain.scale)
post <- posterior_samples(mc.loc, add_chain = T)
mcmc_dens_overlay(post)post %>%
select(lm1, sill, range, scale1, scale2, scale3, scale4) %>%
mcmc_acf()post %>%
select(lm1, sill, range, scale1, scale2, scale3, scale4) %>%
mcmc_trace()post %>%
select(lm1, sill, range, scale1, scale2, scale3, scale4) %>%
mcmc_dens()post %>% select(lm1) %>% mcmc_trace()library(coda)
coda::raftery.diag(posterior_samples(post))>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3740 3746 0.998
> sill 3 4232 3746 1.130
> range 8 9716 3746 2.590
> smooth <NA> <NA> 3746 NA
> scale1 6 8308 3746 2.220
> scale2 3 4095 3746 1.090
> scale3 6 9070 3746 2.420
> scale4 6 8482 3746 2.260
> scale5 8 9888 3746 2.640
> scale6 6 11691 3746 3.120
> scale7 6 7762 3746 2.070
> scale8 6 7760 3746 2.070
> scale9 3 4520 3746 1.210
> scale10 6 8380 3746 2.240
> scale11 6 9414 3746 2.510
> scale12 4 5038 3746 1.340
> scale13 4 7048 3746 1.880
> scale14 3 4128 3746 1.100
> scale15 3 4558 3746 1.220
> scale16 6 8826 3746 2.360
> scale17 4 8092 3746 2.160
> scale18 8 9270 3746 2.470
> scale19 6 8574 3746 2.290
> scale20 9 13014 3746 3.470
> scale21 3 4337 3746 1.160
> scale22 3 4302 3746 1.150
> scale23 4 4635 3746 1.240
> scale24 6 8820 3746 2.350
> scale25 6 8826 3746 2.360
> scale26 4 4752 3746 1.270
> scale27 6 7818 3746 2.090
> scale28 9 12801 3746 3.420
> scale29 6 8190 3746 2.190
> scale30 4 7702 3746 2.060
> scale31 3 4560 3746 1.220
> scale32 6 7708 3746 2.060
> scale33 6 9016 3746 2.410
> scale34 6 8978 3746 2.400
> scale35 12 14139 3746 3.770
> scale36 6 7626 3746 2.040
> scale37 6 8278 3746 2.210
> scale38 6 7260 3746 1.940
> scale39 3 4373 3746 1.170
> scale40 9 13674 3746 3.650
> chain 69075 69075 3746 18.400
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.scale[,c(-4,-9)]),
# as.mcmc(mc2$chain.scale[,c(-4,-9)]),
# as.mcmc(mc3$chain.scale[,c(-4,-9)])))
#
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.scale[,1]),
# as.mcmc(mc2$chain.scale[,1]),
# as.mcmc(mc3$chain.scale[,1])))
BMu1.mcmc<-mcmc.list(as.mcmc(mc1$chain.scale[,c(-4,-9)]),
as.mcmc(mc2$chain.scale[,c(-4,-9)]),
as.mcmc(mc3$chain.scale[,c(-4,-9)]))
summary(BMu1.mcmc)>
> Iterations = 1:10000
> Thinning interval = 1
> Number of chains = 3
> Sample size per chain = 10000
>
> 1. Empirical mean and standard deviation for each variable,
> plus standard error of the mean:
>
> Mean SD Naive SE Time-series SE
> lm1 3.519 1.3910 0.0080307 0.010806
> sill 1.414 1.1966 0.0069083 0.011339
> range 10.441 6.4762 0.0373902 0.063922
> scale1 1.809 0.1598 0.0009226 0.001661
> scale2 2.207 0.1706 0.0009850 0.001530
> scale3 2.088 0.1744 0.0010071 0.001749
> scale4 1.780 0.1572 0.0009075 0.001829
> scale6 2.330 0.2129 0.0012291 0.002264
> scale7 2.901 0.2497 0.0014415 0.002190
> scale8 2.107 0.1737 0.0010027 0.001542
> scale9 2.060 0.1753 0.0010122 0.001552
> scale10 2.136 0.1762 0.0010176 0.001538
> scale11 1.610 0.1552 0.0008963 0.001704
> scale12 2.788 0.3168 0.0018289 0.002945
> scale13 1.932 0.1785 0.0010305 0.001501
> scale14 2.227 0.1890 0.0010912 0.001587
> scale15 1.870 0.1880 0.0010855 0.001704
> scale16 2.003 0.2167 0.0012511 0.001773
> scale17 1.875 0.1783 0.0010294 0.001582
> scale18 2.228 0.2141 0.0012361 0.002856
> scale19 1.689 0.1541 0.0008896 0.001689
> scale20 1.597 0.1923 0.0011104 0.003180
> scale21 2.187 0.2067 0.0011936 0.001578
> scale22 1.791 0.1940 0.0011203 0.001784
> scale23 1.895 0.2030 0.0011722 0.002175
> scale24 2.380 0.2196 0.0012676 0.002243
> scale25 1.549 0.1516 0.0008752 0.001925
> scale26 2.051 0.1780 0.0010275 0.001535
> scale27 1.868 0.1607 0.0009279 0.001863
> scale28 1.538 0.1741 0.0010051 0.002257
> scale29 2.827 0.2687 0.0015512 0.002566
> scale30 2.618 0.2675 0.0015446 0.002603
> scale31 2.083 0.1739 0.0010040 0.001541
> scale32 2.167 0.1487 0.0008585 0.001525
> scale33 1.542 0.1595 0.0009211 0.001994
> scale34 2.125 0.1563 0.0009027 0.001479
> scale35 1.596 0.1905 0.0011001 0.003125
> scale36 1.871 0.1992 0.0011501 0.001985
> scale37 2.101 0.2121 0.0012246 0.002012
> scale38 2.037 0.1885 0.0010886 0.001756
> scale39 1.746 0.1884 0.0010876 0.001579
> scale40 1.526 0.1812 0.0010464 0.002234
> chain 2.000 0.8165 0.0047141 0.000000
>
> 2. Quantiles for each variable:
>
> 2.5% 25% 50% 75% 97.5%
> lm1 1.6039 2.5820 3.227 4.120 7.058
> sill 0.3121 0.6665 1.057 1.711 4.812
> range 2.3281 5.7278 9.019 13.476 26.816
> scale1 1.5099 1.7028 1.806 1.911 2.133
> scale2 1.8930 2.0905 2.198 2.315 2.568
> scale3 1.7694 1.9665 2.081 2.205 2.450
> scale4 1.4845 1.6729 1.778 1.881 2.097
> scale6 1.9246 2.1867 2.323 2.465 2.761
> scale7 2.4643 2.7240 2.887 3.055 3.448
> scale8 1.7847 1.9877 2.104 2.220 2.459
> scale9 1.7218 1.9424 2.059 2.177 2.406
> scale10 1.7938 2.0180 2.132 2.252 2.497
> scale11 1.3121 1.5065 1.608 1.712 1.918
> scale12 2.2050 2.5693 2.769 2.989 3.448
> scale13 1.5807 1.8126 1.933 2.052 2.281
> scale14 1.8770 2.0964 2.220 2.352 2.617
> scale15 1.5099 1.7473 1.865 1.994 2.248
> scale16 1.6073 1.8528 1.996 2.145 2.450
> scale17 1.5136 1.7569 1.878 1.992 2.225
> scale18 1.8129 2.0853 2.223 2.368 2.653
> scale19 1.3903 1.5835 1.687 1.789 1.998
> scale20 1.2398 1.4654 1.593 1.721 2.000
> scale21 1.7864 2.0488 2.185 2.325 2.597
> scale22 1.4287 1.6574 1.785 1.916 2.183
> scale23 1.5234 1.7539 1.883 2.027 2.315
> scale24 1.9740 2.2292 2.371 2.522 2.826
> scale25 1.2591 1.4452 1.545 1.649 1.856
> scale26 1.7034 1.9331 2.047 2.167 2.405
> scale27 1.5682 1.7597 1.863 1.971 2.194
> scale28 1.2143 1.4180 1.531 1.650 1.900
> scale29 2.3530 2.6389 2.807 2.995 3.407
> scale30 2.1277 2.4370 2.604 2.787 3.188
> scale31 1.7489 1.9637 2.082 2.200 2.430
> scale32 1.8911 2.0638 2.164 2.262 2.466
> scale33 1.2348 1.4332 1.538 1.647 1.868
> scale34 1.8289 2.0176 2.122 2.228 2.441
> scale35 1.2396 1.4673 1.590 1.717 1.995
> scale36 1.4885 1.7351 1.870 2.004 2.267
> scale37 1.6931 1.9580 2.098 2.241 2.525
> scale38 1.6741 1.9092 2.035 2.162 2.414
> scale39 1.4048 1.6172 1.737 1.865 2.139
> scale40 1.1760 1.4047 1.523 1.643 1.891
> chain 1.0000 1.0000 2.000 3.000 3.000
xyplot(BMu1.mcmc)densityplot(BMu1.mcmc) #Densidadesplot(BMu1.mcmc)# layout(matrix(1:12, 3,4));
traceplot(BMu1.mcmc)#Grafica de autocorrelacion
autocorr.plot(BMu1.mcmc, auto.layout = TRUE, ask =F)coda::acfplot(BMu1.mcmc) # gelman.diag(BMu1.mcmc)
#
# gelman.plot(BMu1.mcmc)
geweke.diag(BMu1.mcmc)> [[1]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range scale1 scale2 scale3 scale4 scale6
> -0.08032 -0.02250 1.17062 -0.50661 1.31192 0.31638 0.87023 -0.04461
> scale7 scale8 scale9 scale10 scale11 scale12 scale13 scale14
> -0.73503 -1.30912 1.91126 2.23194 0.36890 -0.70658 1.04083 1.01627
> scale15 scale16 scale17 scale18 scale19 scale20 scale21 scale22
> 0.68221 0.73472 1.83715 0.09793 -0.28420 -0.44112 1.22216 -3.26133
> scale23 scale24 scale25 scale26 scale27 scale28 scale29 scale30
> -2.80187 0.04558 -0.13999 -0.41964 1.43725 -1.79599 -1.00720 -0.25150
> scale31 scale32 scale33 scale34 scale35 scale36 scale37 scale38
> -0.66973 1.55224 -0.21793 1.63241 -0.59266 0.97689 0.18939 -0.53277
> scale39 scale40 chain
> 0.69119 -0.37436 NaN
>
>
> [[2]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range scale1 scale2 scale3 scale4 scale6
> -0.08032 -0.02250 1.17062 -0.50661 1.31192 0.31638 0.87023 -0.04461
> scale7 scale8 scale9 scale10 scale11 scale12 scale13 scale14
> -0.73503 -1.30912 1.91126 2.23194 0.36890 -0.70658 1.04083 1.01627
> scale15 scale16 scale17 scale18 scale19 scale20 scale21 scale22
> 0.68221 0.73472 1.83715 0.09793 -0.28420 -0.44112 1.22216 -3.26133
> scale23 scale24 scale25 scale26 scale27 scale28 scale29 scale30
> -2.80187 0.04558 -0.13999 -0.41964 1.43725 -1.79599 -1.00720 -0.25150
> scale31 scale32 scale33 scale34 scale35 scale36 scale37 scale38
> -0.66973 1.55224 -0.21793 1.63241 -0.59266 0.97689 0.18939 -0.53277
> scale39 scale40 chain
> 0.69119 -0.37436 NaN
>
>
> [[3]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range scale1 scale2 scale3 scale4 scale6
> -0.08032 -0.02250 1.17062 -0.50661 1.31192 0.31638 0.87023 -0.04461
> scale7 scale8 scale9 scale10 scale11 scale12 scale13 scale14
> -0.73503 -1.30912 1.91126 2.23194 0.36890 -0.70658 1.04083 1.01627
> scale15 scale16 scale17 scale18 scale19 scale20 scale21 scale22
> 0.68221 0.73472 1.83715 0.09793 -0.28420 -0.44112 1.22216 -3.26133
> scale23 scale24 scale25 scale26 scale27 scale28 scale29 scale30
> -2.80187 0.04558 -0.13999 -0.41964 1.43725 -1.79599 -1.00720 -0.25150
> scale31 scale32 scale33 scale34 scale35 scale36 scale37 scale38
> -0.66973 1.55224 -0.21793 1.63241 -0.59266 0.97689 0.18939 -0.53277
> scale39 scale40 chain
> 0.69119 -0.37436 NaN
#geweke.plot(BMu1.mcmc)
raftery.diag(BMu1.mcmc)> [[1]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3741 3746 0.999
> sill 3 4232 3746 1.130
> range 5 5721 3746 1.530
> scale1 4 4752 3746 1.270
> scale2 3 4095 3746 1.090
> scale3 4 4635 3746 1.240
> scale4 5 5390 3746 1.440
> scale6 4 4674 3746 1.250
> scale7 3 4320 3746 1.150
> scale8 4 4635 3746 1.240
> scale9 3 4520 3746 1.210
> scale10 6 8382 3746 2.240
> scale11 6 9414 3746 2.510
> scale12 4 5038 3746 1.340
> scale13 4 4635 3746 1.240
> scale14 3 4129 3746 1.100
> scale15 3 4558 3746 1.220
> scale16 3 4520 3746 1.210
> scale17 3 4410 3746 1.180
> scale18 8 9272 3746 2.480
> scale19 6 8574 3746 2.290
> scale20 6 9466 3746 2.530
> scale21 3 4338 3746 1.160
> scale22 3 4302 3746 1.150
> scale23 4 4635 3746 1.240
> scale24 4 4635 3746 1.240
> scale25 6 8828 3746 2.360
> scale26 4 4752 3746 1.270
> scale27 6 7818 3746 2.090
> scale28 8 8872 3746 2.370
> scale29 3 4302 3746 1.150
> scale30 3 4483 3746 1.200
> scale31 3 4564 3746 1.220
> scale32 3 4410 3746 1.180
> scale33 6 9018 3746 2.410
> scale34 3 4410 3746 1.180
> scale35 10 12532 3746 3.350
> scale36 4 4832 3746 1.290
> scale37 4 4872 3746 1.300
> scale38 4 4674 3746 1.250
> scale39 3 4374 3746 1.170
> scale40 6 10278 3746 2.740
> chain <NA> <NA> 3746 NA
>
>
> [[2]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3741 3746 0.999
> sill 3 4232 3746 1.130
> range 5 5721 3746 1.530
> scale1 4 4752 3746 1.270
> scale2 3 4095 3746 1.090
> scale3 4 4635 3746 1.240
> scale4 5 5390 3746 1.440
> scale6 4 4674 3746 1.250
> scale7 3 4320 3746 1.150
> scale8 4 4635 3746 1.240
> scale9 3 4520 3746 1.210
> scale10 6 8382 3746 2.240
> scale11 6 9414 3746 2.510
> scale12 4 5038 3746 1.340
> scale13 4 4635 3746 1.240
> scale14 3 4129 3746 1.100
> scale15 3 4558 3746 1.220
> scale16 3 4520 3746 1.210
> scale17 3 4410 3746 1.180
> scale18 8 9272 3746 2.480
> scale19 6 8574 3746 2.290
> scale20 6 9466 3746 2.530
> scale21 3 4338 3746 1.160
> scale22 3 4302 3746 1.150
> scale23 4 4635 3746 1.240
> scale24 4 4635 3746 1.240
> scale25 6 8828 3746 2.360
> scale26 4 4752 3746 1.270
> scale27 6 7818 3746 2.090
> scale28 8 8872 3746 2.370
> scale29 3 4302 3746 1.150
> scale30 3 4483 3746 1.200
> scale31 3 4564 3746 1.220
> scale32 3 4410 3746 1.180
> scale33 6 9018 3746 2.410
> scale34 3 4410 3746 1.180
> scale35 10 12532 3746 3.350
> scale36 4 4832 3746 1.290
> scale37 4 4872 3746 1.300
> scale38 4 4674 3746 1.250
> scale39 3 4374 3746 1.170
> scale40 6 10278 3746 2.740
> chain <NA> <NA> 3746 NA
>
>
> [[3]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 2 3741 3746 0.999
> sill 3 4232 3746 1.130
> range 5 5721 3746 1.530
> scale1 4 4752 3746 1.270
> scale2 3 4095 3746 1.090
> scale3 4 4635 3746 1.240
> scale4 5 5390 3746 1.440
> scale6 4 4674 3746 1.250
> scale7 3 4320 3746 1.150
> scale8 4 4635 3746 1.240
> scale9 3 4520 3746 1.210
> scale10 6 8382 3746 2.240
> scale11 6 9414 3746 2.510
> scale12 4 5038 3746 1.340
> scale13 4 4635 3746 1.240
> scale14 3 4129 3746 1.100
> scale15 3 4558 3746 1.220
> scale16 3 4520 3746 1.210
> scale17 3 4410 3746 1.180
> scale18 8 9272 3746 2.480
> scale19 6 8574 3746 2.290
> scale20 6 9466 3746 2.530
> scale21 3 4338 3746 1.160
> scale22 3 4302 3746 1.150
> scale23 4 4635 3746 1.240
> scale24 4 4635 3746 1.240
> scale25 6 8828 3746 2.360
> scale26 4 4752 3746 1.270
> scale27 6 7818 3746 2.090
> scale28 8 8872 3746 2.370
> scale29 3 4302 3746 1.150
> scale30 3 4483 3746 1.200
> scale31 3 4564 3746 1.220
> scale32 3 4410 3746 1.180
> scale33 6 9018 3746 2.410
> scale34 3 4410 3746 1.180
> scale35 10 12532 3746 3.350
> scale36 4 4832 3746 1.290
> scale37 4 4872 3746 1.300
> scale38 4 4674 3746 1.250
> scale39 3 4374 3746 1.170
> scale40 6 10278 3746 2.740
> chain <NA> <NA> 3746 NA
# heidel.diag(BMu1.mcmc)mc1$chain.shape <- as.data.frame(mc1$chain.shape) %>%
mutate(chain = 1)
mc2$chain.shape <- as.data.frame(mc2$chain.shape) %>%
mutate(chain = 2)
mc3$chain.shape <- as.data.frame(mc3$chain.shape) %>%
mutate(chain = 3)
mc.loc <- mc1$chain.shape %>%
union_all(mc2$chain.shape) %>% union_all(mc3$chain.shape)
post <- posterior_samples(mc.loc, add_chain = T)
mcmc_dens_overlay(post)post %>%
select(lm1, sill, range, shape1, shape2, shape3, shape4) %>%
mcmc_acf()post %>%
select(lm1, sill, range, shape1, shape2, shape3, shape4) %>%
mcmc_trace()post %>%
select(lm1, sill, range, shape1, shape2, shape3, shape4) %>%
mcmc_dens()post %>% select(lm1) %>% mcmc_trace()library(coda)
coda::raftery.diag(posterior_samples(post))>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 4 8092 3746 2.16
> sill 6 8524 3746 2.28
> range 6 8566 3746 2.29
> smooth <NA> <NA> 3746 NA
> shape1 3 4095 3746 1.09
> shape2 3 4095 3746 1.09
> shape3 4 8144 3746 2.17
> shape4 3 4163 3746 1.11
> shape5 3 4520 3746 1.21
> shape6 3 4095 3746 1.09
> shape7 3 4029 3746 1.08
> shape8 6 8532 3746 2.28
> shape9 3 4095 3746 1.09
> shape10 6 8342 3746 2.23
> shape11 4 6654 3746 1.78
> shape12 3 4267 3746 1.14
> shape13 3 4061 3746 1.08
> shape14 3 4267 3746 1.14
> shape15 3 4061 3746 1.08
> shape16 2 3897 3746 1.04
> shape17 2 3961 3746 1.06
> shape18 6 11313 3746 3.02
> shape19 3 4061 3746 1.08
> shape20 4 8092 3746 2.16
> shape21 3 4061 3746 1.08
> shape22 3 4028 3746 1.08
> shape23 3 4095 3746 1.09
> shape24 3 4267 3746 1.14
> shape25 3 4157 3746 1.11
> shape26 2 3865 3746 1.03
> shape27 2 3929 3746 1.05
> shape28 6 7792 3746 2.08
> shape29 3 4095 3746 1.09
> shape30 2 3994 3746 1.07
> shape31 3 4410 3746 1.18
> shape32 3 4028 3746 1.08
> shape33 3 4095 3746 1.09
> shape34 4 8298 3746 2.22
> shape35 6 8736 3746 2.33
> shape36 2 3994 3746 1.07
> shape37 3 4128 3746 1.10
> shape38 3 4128 3746 1.10
> shape39 3 4061 3746 1.08
> shape40 3 4095 3746 1.09
> chain 69075 69075 3746 18.40
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.shape[,c(-4,-9)]),
# as.mcmc(mc2$chain.shape[,c(-4,-9)]),
# as.mcmc(mc3$chain.shape[,c(-4,-9)])))
#
# coda::gelman.diag(mcmc.list(as.mcmc(mc1$chain.shape[,1]),
# as.mcmc(mc2$chain.shape[,1]),
# as.mcmc(mc3$chain.shape[,1])))
BMu1.mcmc<-mcmc.list(as.mcmc(mc1$chain.shape[,c(-4,-9)]),
as.mcmc(mc2$chain.shape[,c(-4,-9)]),
as.mcmc(mc3$chain.shape[,c(-4,-9)]))
summary(BMu1.mcmc)>
> Iterations = 1:10000
> Thinning interval = 1
> Number of chains = 3
> Sample size per chain = 10000
>
> 1. Empirical mean and standard deviation for each variable,
> plus standard error of the mean:
>
> Mean SD Naive SE Time-series SE
> lm1 0.7418 0.42890 0.0024763 0.0045351
> sill 0.1937 0.32464 0.0018743 0.0050592
> range 3.5989 5.38595 0.0310958 0.0908870
> shape1 0.6459 0.09745 0.0005626 0.0007165
> shape2 0.2402 0.10805 0.0006238 0.0009879
> shape3 0.6073 0.09190 0.0005306 0.0006901
> shape4 0.6922 0.09533 0.0005504 0.0007864
> shape6 0.6417 0.10743 0.0006202 0.0008532
> shape7 0.4316 0.08297 0.0004790 0.0005468
> shape8 0.5020 0.09171 0.0005295 0.0006737
> shape9 0.4534 0.10082 0.0005821 0.0007489
> shape10 0.5147 0.10265 0.0005926 0.0007973
> shape11 0.6614 0.09924 0.0005730 0.0007437
> shape12 0.8658 0.13977 0.0008070 0.0012274
> shape13 0.4066 0.10210 0.0005895 0.0007181
> shape14 0.4597 0.10748 0.0006205 0.0007699
> shape15 0.5024 0.10674 0.0006163 0.0007605
> shape16 0.4802 0.10774 0.0006221 0.0007219
> shape17 0.5337 0.11206 0.0006470 0.0007879
> shape18 0.6569 0.11091 0.0006403 0.0008731
> shape19 0.6658 0.09729 0.0005617 0.0007213
> shape20 0.9772 0.11750 0.0006784 0.0012123
> shape21 0.5361 0.10946 0.0006320 0.0008032
> shape22 0.6036 0.11916 0.0006880 0.0009378
> shape23 0.5628 0.12417 0.0007169 0.0010802
> shape24 0.5525 0.10463 0.0006041 0.0007478
> shape25 0.7334 0.10058 0.0005807 0.0008044
> shape26 0.5391 0.10400 0.0006004 0.0007936
> shape27 0.6670 0.09224 0.0005325 0.0007783
> shape28 0.8440 0.10508 0.0006067 0.0009787
> shape29 0.3711 0.10235 0.0005909 0.0007183
> shape30 0.5511 0.11441 0.0006605 0.0008423
> shape31 0.5225 0.10638 0.0006142 0.0008025
> shape32 0.4321 0.06945 0.0004010 0.0005309
> shape33 0.7575 0.09803 0.0005660 0.0007805
> shape34 0.4100 0.08273 0.0004776 0.0005993
> shape35 0.9596 0.11721 0.0006767 0.0011816
> shape36 0.6524 0.11323 0.0006537 0.0008618
> shape37 0.6166 0.11355 0.0006556 0.0008263
> shape38 0.5961 0.09975 0.0005759 0.0007665
> shape39 0.3269 0.11645 0.0006723 0.0009777
> shape40 0.8475 0.11604 0.0006700 0.0009448
> chain 2.0000 0.81651 0.0047141 0.0000000
>
> 2. Quantiles for each variable:
>
> 2.5% 25% 50% 75% 97.5%
> lm1 0.21272 0.54460 0.64613 0.8139 1.8975
> sill 0.02498 0.05159 0.08874 0.1944 1.0693
> range 0.25977 0.76101 1.64603 4.0021 19.5252
> shape1 0.46191 0.57955 0.64278 0.7097 0.8461
> shape2 0.02828 0.16409 0.24133 0.3159 0.4457
> shape3 0.43191 0.54578 0.60458 0.6658 0.7940
> shape4 0.51797 0.62601 0.68734 0.7506 0.8940
> shape6 0.44422 0.56619 0.63709 0.7135 0.8661
> shape7 0.27640 0.37399 0.42936 0.4868 0.6012
> shape8 0.32423 0.44009 0.50070 0.5637 0.6837
> shape9 0.25684 0.38609 0.45436 0.5197 0.6524
> shape10 0.32920 0.44343 0.50876 0.5784 0.7334
> shape11 0.46621 0.59598 0.66084 0.7268 0.8572
> shape12 0.61774 0.76700 0.85804 0.9559 1.1608
> shape13 0.21081 0.33785 0.40406 0.4735 0.6119
> shape14 0.25519 0.38571 0.46042 0.5305 0.6767
> shape15 0.30317 0.42997 0.49838 0.5696 0.7258
> shape16 0.28151 0.40469 0.47570 0.5500 0.7029
> shape17 0.32500 0.45801 0.53085 0.6045 0.7658
> shape18 0.45875 0.57873 0.65075 0.7275 0.8896
> shape19 0.47342 0.60086 0.66421 0.7317 0.8617
> shape20 0.76663 0.89403 0.97158 1.0527 1.2233
> shape21 0.33952 0.45975 0.52967 0.6072 0.7665
> shape22 0.37798 0.52251 0.60138 0.6814 0.8406
> shape23 0.32039 0.47840 0.56323 0.6486 0.8061
> shape24 0.34594 0.48262 0.55170 0.6218 0.7588
> shape25 0.53804 0.66647 0.73268 0.7997 0.9355
> shape26 0.34906 0.46804 0.53500 0.6029 0.7594
> shape27 0.49903 0.60412 0.66181 0.7250 0.8598
> shape28 0.64710 0.77226 0.84098 0.9125 1.0580
> shape29 0.17279 0.30149 0.37058 0.4387 0.5759
> shape30 0.33797 0.47217 0.54881 0.6280 0.7807
> shape31 0.31713 0.45309 0.52170 0.5908 0.7373
> shape32 0.29837 0.38406 0.43187 0.4788 0.5692
> shape33 0.57353 0.69037 0.75296 0.8216 0.9599
> shape34 0.25091 0.35412 0.40783 0.4646 0.5771
> shape35 0.74463 0.87792 0.95493 1.0362 1.2030
> shape36 0.44284 0.57410 0.64716 0.7265 0.8855
> shape37 0.40489 0.53909 0.61311 0.6907 0.8530
> shape38 0.41228 0.52667 0.59225 0.6608 0.8029
> shape39 0.10069 0.24825 0.32669 0.4045 0.5541
> shape40 0.63246 0.76711 0.84226 0.9239 1.0875
> chain 1.00000 1.00000 2.00000 3.0000 3.0000
xyplot(BMu1.mcmc)densityplot(BMu1.mcmc) #Densidadesplot(BMu1.mcmc)# layout(matrix(1:12, 3,4));
traceplot(BMu1.mcmc)#Grafica de autocorrelacion
autocorr.plot(BMu1.mcmc, auto.layout = TRUE, ask =F)coda::acfplot(BMu1.mcmc) # gelman.diag(BMu1.mcmc)
#
# gelman.plot(BMu1.mcmc)
geweke.diag(BMu1.mcmc)> [[1]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range shape1 shape2 shape3 shape4 shape6
> -0.541921 -0.911255 -0.635529 0.546572 -0.075900 -0.584468 1.271316 -0.874372
> shape7 shape8 shape9 shape10 shape11 shape12 shape13 shape14
> 1.030963 -0.971523 -0.935714 0.256445 2.004252 0.695535 -0.413436 0.038511
> shape15 shape16 shape17 shape18 shape19 shape20 shape21 shape22
> 1.790014 0.875726 0.739533 0.005974 1.728120 0.420641 -0.812484 -0.303818
> shape23 shape24 shape25 shape26 shape27 shape28 shape29 shape30
> -1.669508 -0.062697 1.212304 -0.195292 0.512424 1.141296 0.221783 1.097566
> shape31 shape32 shape33 shape34 shape35 shape36 shape37 shape38
> 0.015901 0.606067 2.901314 0.798967 0.277861 -0.220946 -0.381475 0.982637
> shape39 shape40 chain
> 0.745817 0.855370 NaN
>
>
> [[2]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range shape1 shape2 shape3 shape4 shape6
> -0.541921 -0.911255 -0.635529 0.546572 -0.075900 -0.584468 1.271316 -0.874372
> shape7 shape8 shape9 shape10 shape11 shape12 shape13 shape14
> 1.030963 -0.971523 -0.935714 0.256445 2.004252 0.695535 -0.413436 0.038511
> shape15 shape16 shape17 shape18 shape19 shape20 shape21 shape22
> 1.790014 0.875726 0.739533 0.005974 1.728120 0.420641 -0.812484 -0.303818
> shape23 shape24 shape25 shape26 shape27 shape28 shape29 shape30
> -1.669508 -0.062697 1.212304 -0.195292 0.512424 1.141296 0.221783 1.097566
> shape31 shape32 shape33 shape34 shape35 shape36 shape37 shape38
> 0.015901 0.606067 2.901314 0.798967 0.277861 -0.220946 -0.381475 0.982637
> shape39 shape40 chain
> 0.745817 0.855370 NaN
>
>
> [[3]]
>
> Fraction in 1st window = 0.1
> Fraction in 2nd window = 0.5
>
> lm1 sill range shape1 shape2 shape3 shape4 shape6
> -0.541921 -0.911255 -0.635529 0.546572 -0.075900 -0.584468 1.271316 -0.874372
> shape7 shape8 shape9 shape10 shape11 shape12 shape13 shape14
> 1.030963 -0.971523 -0.935714 0.256445 2.004252 0.695535 -0.413436 0.038511
> shape15 shape16 shape17 shape18 shape19 shape20 shape21 shape22
> 1.790014 0.875726 0.739533 0.005974 1.728120 0.420641 -0.812484 -0.303818
> shape23 shape24 shape25 shape26 shape27 shape28 shape29 shape30
> -1.669508 -0.062697 1.212304 -0.195292 0.512424 1.141296 0.221783 1.097566
> shape31 shape32 shape33 shape34 shape35 shape36 shape37 shape38
> 0.015901 0.606067 2.901314 0.798967 0.277861 -0.220946 -0.381475 0.982637
> shape39 shape40 chain
> 0.745817 0.855370 NaN
#geweke.plot(BMu1.mcmc)
raftery.diag(BMu1.mcmc)> [[1]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 3 4163 3746 1.11
> sill 4 5080 3746 1.36
> range 5 6077 3746 1.62
> shape1 3 4095 3746 1.09
> shape2 3 4095 3746 1.09
> shape3 3 4338 3746 1.16
> shape4 3 4163 3746 1.11
> shape6 3 4095 3746 1.09
> shape7 3 4031 3746 1.08
> shape8 2 3865 3746 1.03
> shape9 3 4095 3746 1.09
> shape10 2 3962 3746 1.06
> shape11 2 3962 3746 1.06
> shape12 3 4267 3746 1.14
> shape13 3 4061 3746 1.08
> shape14 3 4267 3746 1.14
> shape15 3 4061 3746 1.08
> shape16 2 3897 3746 1.04
> shape17 2 3962 3746 1.06
> shape18 2 3929 3746 1.05
> shape19 3 4061 3746 1.08
> shape20 4 4752 3746 1.27
> shape21 3 4061 3746 1.08
> shape22 3 4028 3746 1.08
> shape23 3 4095 3746 1.09
> shape24 3 4267 3746 1.14
> shape25 3 4147 3746 1.11
> shape26 2 3865 3746 1.03
> shape27 2 3929 3746 1.05
> shape28 3 4113 3746 1.10
> shape29 3 4095 3746 1.09
> shape30 2 3994 3746 1.07
> shape31 3 4410 3746 1.18
> shape32 3 4028 3746 1.08
> shape33 3 4095 3746 1.09
> shape34 3 4129 3746 1.10
> shape35 6 8738 3746 2.33
> shape36 2 3994 3746 1.07
> shape37 3 4129 3746 1.10
> shape38 3 4129 3746 1.10
> shape39 3 4061 3746 1.08
> shape40 3 4095 3746 1.09
> chain <NA> <NA> 3746 NA
>
>
> [[2]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 3 4163 3746 1.11
> sill 4 5080 3746 1.36
> range 5 6077 3746 1.62
> shape1 3 4095 3746 1.09
> shape2 3 4095 3746 1.09
> shape3 3 4338 3746 1.16
> shape4 3 4163 3746 1.11
> shape6 3 4095 3746 1.09
> shape7 3 4031 3746 1.08
> shape8 2 3865 3746 1.03
> shape9 3 4095 3746 1.09
> shape10 2 3962 3746 1.06
> shape11 2 3962 3746 1.06
> shape12 3 4267 3746 1.14
> shape13 3 4061 3746 1.08
> shape14 3 4267 3746 1.14
> shape15 3 4061 3746 1.08
> shape16 2 3897 3746 1.04
> shape17 2 3962 3746 1.06
> shape18 2 3929 3746 1.05
> shape19 3 4061 3746 1.08
> shape20 4 4752 3746 1.27
> shape21 3 4061 3746 1.08
> shape22 3 4028 3746 1.08
> shape23 3 4095 3746 1.09
> shape24 3 4267 3746 1.14
> shape25 3 4147 3746 1.11
> shape26 2 3865 3746 1.03
> shape27 2 3929 3746 1.05
> shape28 3 4113 3746 1.10
> shape29 3 4095 3746 1.09
> shape30 2 3994 3746 1.07
> shape31 3 4410 3746 1.18
> shape32 3 4028 3746 1.08
> shape33 3 4095 3746 1.09
> shape34 3 4129 3746 1.10
> shape35 6 8738 3746 2.33
> shape36 2 3994 3746 1.07
> shape37 3 4129 3746 1.10
> shape38 3 4129 3746 1.10
> shape39 3 4061 3746 1.08
> shape40 3 4095 3746 1.09
> chain <NA> <NA> 3746 NA
>
>
> [[3]]
>
> Quantile (q) = 0.025
> Accuracy (r) = +/- 0.005
> Probability (s) = 0.95
>
> Burn-in Total Lower bound Dependence
> (M) (N) (Nmin) factor (I)
> lm1 3 4163 3746 1.11
> sill 4 5080 3746 1.36
> range 5 6077 3746 1.62
> shape1 3 4095 3746 1.09
> shape2 3 4095 3746 1.09
> shape3 3 4338 3746 1.16
> shape4 3 4163 3746 1.11
> shape6 3 4095 3746 1.09
> shape7 3 4031 3746 1.08
> shape8 2 3865 3746 1.03
> shape9 3 4095 3746 1.09
> shape10 2 3962 3746 1.06
> shape11 2 3962 3746 1.06
> shape12 3 4267 3746 1.14
> shape13 3 4061 3746 1.08
> shape14 3 4267 3746 1.14
> shape15 3 4061 3746 1.08
> shape16 2 3897 3746 1.04
> shape17 2 3962 3746 1.06
> shape18 2 3929 3746 1.05
> shape19 3 4061 3746 1.08
> shape20 4 4752 3746 1.27
> shape21 3 4061 3746 1.08
> shape22 3 4028 3746 1.08
> shape23 3 4095 3746 1.09
> shape24 3 4267 3746 1.14
> shape25 3 4147 3746 1.11
> shape26 2 3865 3746 1.03
> shape27 2 3929 3746 1.05
> shape28 3 4113 3746 1.10
> shape29 3 4095 3746 1.09
> shape30 2 3994 3746 1.07
> shape31 3 4410 3746 1.18
> shape32 3 4028 3746 1.08
> shape33 3 4095 3746 1.09
> shape34 3 4129 3746 1.10
> shape35 6 8738 3746 2.33
> shape36 2 3994 3746 1.07
> shape37 3 4129 3746 1.10
> shape38 3 4129 3746 1.10
> shape39 3 4061 3746 1.08
> shape40 3 4095 3746 1.09
> chain <NA> <NA> 3746 NA
# heidel.diag(BMu1.mcmc)bg <- bg %>%
dplyr::filter(Rain > 30)
bg %>%
plot_density(x = Rain,
title = "(2018-2020) basic density plot general of cum_rain",
fill = "salmon1")
table(bg$Station)
library(SpatialExtremes)
## Not run:
## Generate realizations from the model
n.site <- 4 # numero de sitios muestreados
n.obs <- 11 # numero de observaciones tomada en cada sitio muestreado
coord <- bg %>% # creamos la matrix de coordenadas
select(Longitud, Latitud) %>% # esta matrix de coordenadas solo incluye el
unique() # numero de sitios donde se muestreo
coord <- cbind(lon = coord$Longitud, lat = coord$Latitud) # matrix de coordenadas
# sin mediciones repetidas
# class(coord)
# coord <- cbind(lon = runif(n.site, -10, 10), lat = runif(n.site, -10 , 10))
# creamos un proceso gaussiano para cada parametro de la GEV, esta parte se deja
# a elección del autor elegirlo.
gp.loc <- SpatialExtremes::rgp(1, coord, "powexp", sill = 0.01, range = 10, smooth = 1)
gp.scale <- SpatialExtremes::rgp(1, coord, "powexp", sill = 0.001, range = 15, smooth = 1)
gp.shape <- SpatialExtremes::rgp(1, coord, "powexp", sill = 0.001, range = 20, smooth = 1)
locs <- 26 + 0.5 * coord[,"lon"] + gp.loc
scales <- 10 + 0.2 * coord[,"lat"] + gp.scale
shapes <- 0.15 + gp.shape
SpatialExtremes::rgev(n.obs, locs[1], scales[1], shapes[1])
# Estimated Parameters:
# xi mu beta
# 0.3944641 34.1387472 4.7882714
data <- matrix(NA, n.obs, n.site)
#
# for (i in 1:n.site)
# data[,i] <- SpatialExtremes::rgev(n.obs, locs[i], scales[i], shapes[i])
dim(bg)
table(bg$Station)
data[,1] <-bg %>%
dplyr::select(Rain, Station) %>%
dplyr::filter(Station == "Arriaga") %>%
# slice(1:5) %>%
dplyr::select(Rain) %>%
dplyr::arrange(-Rain) %>%
slice(1:11) %>%
as.matrix()
data[,2] <-bg %>%
dplyr::select(Rain, Station) %>%
dplyr::filter(Station == "Comitan") %>%
# slice(1:5) %>%
select(Rain) %>%
arrange(-Rain) %>% slice(1:11) %>%
as.matrix()
# data[,3] <- bg %>%
# dplyr::select(cum_rain, Station) %>%
# dplyr::filter(Station == "scdlc") %>%
# # slice(1:5) %>%
# select(cum_rain) %>%
# arrange(-cum_rain) %>% # slice(1:61) %>%
# as.matrix()
data[,3] <- bg %>%
dplyr::select(Rain, Station) %>%
dplyr::filter(Station == "Tapachula") %>%
# slice(1:5) %>%
select(Rain) %>%
arrange(-Rain) %>% slice(1:11) %>%
as.matrix()
data[,4] <- bg %>%
dplyr::select(Rain, Station) %>%
dplyr::filter(Station == "Tuxtla") %>%
# slice(1:5) %>%
select(Rain) %>%
arrange(-Rain) %>% slice(1:11) %>%
as.matrix()
# b2018[b2018$Station == "Arriaga", 3]
# data[,5] <- b2018[b2018$Station == "tuxtlag", 3]
# data[,1] <- b2018[1:10,3]
# loc.form <- y ~ lon
# scale.form <- y ~ lat
# shape.form <- y ~ 1
loc.form <- y ~ 1
scale.form <- y ~ 1
shape.form <- y ~ 1
hyper <- list()
hyper$sills <- list(loc = c(1,8),
scale = c(1,1),
shape = c(1,0.02))
hyper$ranges <- list(loc = c(2,20),
scale = c(1,5),
shape = c(1, 10))
hyper$smooths <- list(loc = c(1,1/3),
scale = c(1,1/3),
shape = c(1, 1/3))
hyper$betaMeans <- list(loc = 9,
scale = 6,
shape = 2)
hyper$betaIcov <- list(loc = solve(diag(c(10), 1, 1)),
scale = solve(diag(c(10), 1, 1)),
shape = solve(diag(c(0.13), 1, 1)))
## We will use an exponential covariance function so the jump sizes for
## the shape parameter of the covariance function are null.
prop <- list(gev = c(1.2, 0.08, 0.08), ranges = c(0.7, 0.8, 0.7), smooths = c(0,0,0))
start <- list(sills = c(4, .36, 0.009),
ranges = c(24, 17, 16),
smooths= c(1, 1, 1),
beta = list(loc = c(24),
scale = c(7.31),
shape = c(0.54)))
# Estimated Parameters:
# xi mu beta
# 0.5304988 37.3221090 7.4043417
mc1 <- latent(data, coord,
loc.form = loc.form,
scale.form = scale.form,
shape.form = shape.form,
hyper = hyper,
prop = prop,
start = start,
n = 10000,
burn.in = 5000,
thin = 15)
mc1
mc2 <- latent(data, coord,
loc.form = loc.form,
scale.form = scale.form,
shape.form = shape.form,
hyper = hyper,
prop = prop,
start = start,
n = 10000,
burn.in = 5000,
thin = 15)
mc2
mc3 <- latent(data, coord,
loc.form = loc.form,
scale.form = scale.form,
shape.form = shape.form,
hyper = hyper,
prop = prop,
start = start,
n = 10000,
burn.in = 5000,
thin = 15)
mc3
## End(Not run)
# a <- summary(mc)
#
# a
# head(mc$chain.loc)
# head(mc$chain.scale)
# head(mc$chain.shape)
# head(mc$loc.dsgn.mat)
# head(mc$scale.dsgn.mat)
# head(mc$shape.dsgn.mat)
# tidybayes::add_residual_draws(mc$chain.loc[,6:10])
# mod1_sim <- coda::coda.samples(model = mc,
# variable.names = chain.loc,
# n.iter = 5000)